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We present a coupled atomistic-continuum method for the modeling of defects and interface dynamics 
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conditions between the atomistic and continuum regions. These conditions ensure the accurate passage of large 
scale information between the atomistic and continuum regions and at the same time minimize the reflection 
of phonons at the atomistic-continuum interface. They can be made adaptive if we choose appropriate weight 
functions. We present applications to dislocation dynamics, friction between two-dimensional crystal surfaces 
and fracture dynamics. We compare results of the coupled method and the detailed atomistic model. 
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1. INTRODUCTION 

Traditionally two apparently separate approaches have been used to model a continuous medium. 
The first is the continuum theory, in the form of partial differential equations describing the conser- 
vation laws and constitutive relations. This approach has been impressively successful in a number 
of areas such as solid and fluid mechanics. It is very efficient, simple and often involves very few ma- 
terial parameters. But it becomes inaccurate for problems in which the detailed atomistic processes 
affect the macroscopic behavior of the medium, or when the scale of the medium is small enough 
that the continuum approximation becomes questionable. Such situations are often found in studies 
of properties and defects of micro- or nano- systems and devices. The second approach is atomistic, 
aiming at finding the detailed behavior of each individual atom using molecular dynamics or quantum 
mechanics. This approach can in principle accurately model the underlying physical processes. But 
it is often times prohibitively expensive. 

Recently an alternative approach has been explored that couples the atomistic and continuum 
approaches j^, ^ |^, ^ ^, ^, The main idea is to use atomistic modeling at places where the 
displacement field varies on an atomic scale, and the continuum approach elsewhere. The most 
successful and best-known implementation is the quasi-continuum method |^ which combines the 
adaptive finite element procedure with an atomistic evaluation of the potential energy of the system. 
This method has been applied to a number of examples ^ |l^ , and interesting details were learned 
about the structure of crystal defects. 

Extension of the quasi-continuum method to dynamic problems has not been straightforward 

1^, 0. The main difficulty lies in the proper matching between the atomistic and continuum regions. 
Since the details of lattice vibrations, the phonons, which are an intrinsic part of the atomistic model, 
cannot be represented at the continuum level, conditions must be met that the phonons are not 
reflected at the atomistic-continuum interface. Since the atomistic region is expected to be a very 
small part of the computational domain, violation of this condition quickly leads to local heating of 
the atomistic region and destroys the simulation. In addition, the matching between the atomistic- 
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continuum interface has to be such that large scale information is accurately transmitted in both 
directions. 

The main purpose of the present paper is to introduce a new class of matching conditions between 
atomistic and continuum regions. These matching conditions have the property that they allow 
accurate passage of large scale (scales that are represented by the continuum model) information 
between the atomistic and continuum regions and no reflection of phonon energy to the atomistic 
region. These conditions can also be used in pure molecular dynamics simulations as the border 
conditions to ensure no reflection of phonons at the boundary of the simulation. As applications, 
we use our method to study the dynamics of dislocations in the Frenkel-Kontorova model, friction 
between crystal surfaces and crack propagation. 

2. CONTINUUM APPROXIMATION OF ATOMISTIC MODELS 

As a first step toward constructing a coupled atomistic-continuum method, we discuss briefly how 
continuum equations are obtained from atomistic models. 

2.1. ID Prenkel-Kontorova Model — the Klein-Gordan Equation 

We first consider a simple problem, the Frenkel-Kontorova Model. This is a one-dimensional chain 
of particles in a periodic potential, coupled by springs. We will take the potential to be: 



Here a is the equilibrium distance between neighboring particles, int(x/a) is the integer part oi x/a. 
Denote by Xn the position of the n— th particle, the dynamic equation for the particles is given by 



where / is the applied force. 

One interesting aspect of the Prenkel-Kontorova model is the possibility of having dislocations 
in the system, which corresponds to vacant or doubly occupied potential wells. In the absence of 
dislocations, the equilibrium positions of the particles are given by xj = ja. In general, we let 




(2.1) 



mXn = k[Xn+i -Xn-a\+ k[Xn-l - 



Xn + a]- U'{Xn) + f. 



(2.2) 
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Xj = a{j + Uj) and / = f /a. u is then the displacement field. A dislocation corresponds to a kink in 
u. Far from the dislocations, we can assume \uj — [uj\ \ <^ 1 where [u] is the integer part of u. Then 
we get, assuming [uj] — 0, 

milj = fc[uj_).i — 2uj + Uj-i] — K-Uj + f. (2-3) 
Let r = tyjm/{ka?), K, ~ K,/{ka?), and f ~ f / {ka^), we obtain: 

= Uj+i - 2uj + Uj-i - ICuj + f. (2.4) 
Taking the limit as a ^ 0, we obtain the continuum limit equation for the displacement field u, 

This is simply the Klein-Gordan equation. 

2.2. 2D Triangular Lattice — Isotropic Elasticity 

Now we consider the triangular lattice model. We assume that nearest neighbor atoms interact 
via central forces whose potential is given by $(r'^) where r is the distance between the atoms (see 
Figure [|) . From Newton's law, we have 

mro = -E^ro,<i>(|rojn, (2.6) 
j 

where m is the mass of the atoms, rj is the position of the j-th atom (j = (ji, j^)), fo.j — ro ^ fj- 
Let {Rj} be the equilibrium positions of the atoms. The lattice constant a satisfies the equilibrium 
condition <i>'(a^) = 0. Let {uj} be the displacement vectors, uj = r j — Rj. Taylor expanding and 
omitting nonlinear terms in u, we get 

= -2^<I>'(|rojnroj 
j 

= -2^[$'(a2) + 2a>"(a2)Rj • (uj - Uo)]hRj + Uo - uj] 
j 

= 4ci>"(a2)^(Rj®Rj)(uj-Uo). (2.7) 



Take the example of a Lennard- Jones potential: 



$(r) = eo 



(2.8) 



,(r/ao)i2 (r/ao)6, 

then the lattice constant is equal to under the assumption of nearest neighbor interaction. 



In this case, equation (2/7) becomes 



m.. 18 
— Uo = 



1 

4 



/ 



(ui,o - 2uo,o + u_i,o) 
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%/3 


4 


4 




1 
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4 y 



(uo,i - 2uo,o + Uo -i) 



4 



V 4 4 y 



(U_ij - 2Uo,o + Ui _i) 



(2.9) 



Let T = t^/m/eQ. The continuum limit (as a ^ 0) of the above equations is 
^ ^ ( 27/4 ^ ( 



81/4 

27/4 
(A + Ai)V(V • u) + ^ A u, 



V 



81/4 



V 



54/4 
54/4 



dxdy 



(2.10) 



where X — fi = This is the equation for isotropic elasticity. 



2.3. Slepyan Model of Fracture 

Here we give the one-dimensional and two-dimensional Slepyan models of fracture [l^ . In the ID 
case, one can view it as a model for the atoms lying along a crack surface. Nearest neighbors are 
connected by elastic springs, with spring constant fc, and the atoms are tied to the other side of the 
crack surface by similar springs, which however snap when extended past some breaking point. The 
lines of atoms are being pulled apart by weak springs of spring constant k/N. These weak springs 
are meant schematically to represent N vertical rows of atoms pulling in series. Let {uj^+,Uj^-} be 
the displacement of atoms on the top and bottom crack surfaces respectively. The equation which 
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describes the upper line of mass points in this model is 

fc(iij+i,+ - 2wj,+ + 

muj^+ = < (2.11) 

Here, the first term at the right hand side is elastic coupling to neighbors, the second term is the 
driving force by displacing edges of the strip, the third term is the bonding to atoms at the opposite 
side of the crack surface, the last term is the dissipation, is a step function, and the term containing 
it describes bonds which snap when their total extension reaches a distance 2m/, where is a fracture 
distance. Assume the lattice constant is a. In the region far away from the fracture, we have 

k 

Tnuj,+ = — 2uj^+ + + j;^{Un — + k{uj- — Uj,+) — buj,+ (2-12) 

Dividing by and taking the limit as a — > 0, we obtain 



with T = t^m/{ka?), b = b/ (ka). 

Now we consider a simple 2D model (see Figure ||). A crack moves in a lattice strip composed of 
2A^ rows of mass points. Assume that all the atoms are located at square lattice points if there is no 
exterior force on them. All of the bonds between lattice points are brittle-elastic, behaving as perfect 
linear springs until the instant they snap, from which point they exert no force. The displacement 
of each mass point is described by a single spatial coordinate u^.j, which can be interpreted as the 
height of mass point into or out of the page. The index i takes integer values, while j = 

1/2- N, - ■■ ,-1/2,1/2, - ■ ■,N -1/2. The model is described by the equation 

mili^j ^ -biiij + ^ T{ui>j'-Ui^j), (2.14) 

nearest 
neighbors 

with 

T{r) = kre{2uf - \r\) (2.15) 



representing the brittle nature of the springs, 9 the step function, and h the coefficient of a small 
dissipative term. The boundary condition which drives the motion of the crack is 



Mi,±(Ar-l/2) 



= ±U 



N- 



(2.16) 



Similarly, we can get the continuum limit of (2.14) 



d'^u ^ ~du 
Ot 



(2.17) 



3. PHONONS 



Among the most essential differences between the atomistic and continuum behavior is the presence 
of phonons, the lattice vibrations, at the atomistic scale. In this section we will briefly review the 
spectrum of the phonons. Let us first consider the simplest model: ID discrete wave equation ( ^.3| ) 
with k — 1, K, = Q and f — 0- After discretization in time, we have 

m"+1-2u" + <-i 



At2 



2.- + K" . 



(3.1) 



where is the displacement of the j-th particle at time t = nAt. 

The phonon spectrum for (3.1) is obtained by looking for solutions of the type u" = e*(""'^*+^^' . 
This gives us the dispersion relation 



-— sm — ■ 
At 2 



ujAt . ^ 

— sm — . 

2 



(3.2) 



For the case when At = 0.01, this dispersion relation is depicted in Figure y. If /C 7^ 0, we have 

1 



LjAt 

— — sm — — 
At 2 



sin2|+/C/4. 



(3.3) 



Consider now the 2D triangular lattice described by (2.7). Let us look for the solutions of the type 



= e*(? '"-i"""'^*)U with ^ = (6,6)^- Substituting this expression into (|j|), we obtain 

2 

I U 



sm 



At/2 



18 

AU. 



4 sin 



2 Cia 



,2 + 



a + sm 



•2 Ci-V36 



a V3( 



sm 



2 ;i+V3;2 



a — sm 



2 $i-\/36 . 



V3(sin2 ii±VMia - sin^ «i^^a) 



3 (sin 



2 Ci+\/36 . 



^2 ;i~\/3C2 



a) 



U 



(3.4) 



The eigenvalues of the matrix A are given by, 

A±-^{a + /3 + 7± + /?2 + 72 - a/3 - a7 - /37} , (3.5) 

where 

• 2 6 a -26 + V3^2 • 2 - 

a = sm —a, p = sm a, 7 = sm ; a. 

2 4 4 

The dispersion relation now has two branches 

2 2 2 2 

'^p(Ci,6) = ^arcsin( — v/A+), ^^(^1, 6) = ^ arcsin( — VA-), (3.6) 

where "p" and "s" stands for "pressure" and "shear" waves respectively. 

4. OPTIMAL LOCAL MATCHING CONDITIONS 

We now come to the interface between the atomistic and continuum regions. As we mentioned 
earlier, designing proper matching conditions at this interface is a major challenge in such a coupled 
atomistic/continuum approach. The basic requirements for the matching conditions are the following: 

(1) . Reflection of phonons to the atomistic region should be minimal. This is particularly crucial 
since the atomistic regions are typically very small for the purpose of computational efficiency, reflec- 
tion of phonon energy back to the atomistic region will trigger local heating and melt the crystalline 
structure. 

(2) . Accurate exchange of large scale information between the atomistic and continuum regions. 

The first requirement is reminiscent of the absorbing boundary conditions required for the compu- 
tation of waves H, 0. Indeed our work draws much inspiration from that literature. There are some 
crucial differences between the phonon problem considered here and the ones studied in the litera- 
ture on absorbing boundary conditions. The most obvious one is the fact that the electromagnetic or 
acoustic waves are continuum objects modeled by partial differential equations, and the associated ab- 
sorbing boundary conditions often use small wavenumber and/or frequency approximations, whereas 
the phonons are intrinsically discrete with substantial energy distributed at high wavenumbers. 

In the following we will give an example of a simple discrete wave equation for which exact reflec- 
tionless boundary conditions can be found. Such exact boundary conditions are highly nonlocal and 



therefore not practical. But they give us guidelines on how approximate boundary conditions should 
be constructed. We then present a method that constructs optimal local matching conditions, given 
a predetermined stencil. 



4.1. Exact Boundary Conditions for ID Discrete Wave Equation 



Consider equation (3.1). It is supposed to be solved for all integer values of j. Now let us assume 
that we will truncate the computational domain and only compute u" for j > 0. Assuming there are 
no sources of waves coming from j < 0, we still want to obtain the same solution as if the computation 
is done for all j. At j = 0, we will impose a new boundary condition to make sure that the phonons 
arriving from j > are not reflected back at j = 0. 

At j = 0, we replace (^) by 

< = E ''kju]-\ ao.o = 0. (4.1) 
We would like to determine the coefficients {ak.j}- For the simple problem at hand, it is possible to 



obtain analytical formulas of {akj} such that the imposition of (O) together with the solution of 



(3.1) for j > reproduces exactly the solution of (3J) if it was solved for all integer values of j, i.e 



an exact reflectionless boundary condition can be found. 

First, let us consider the case of /C = and / = 0. Let A = At and let us look for solutions of the 
form: 

< = ^"e, 1^1 < 1- (4.2) 



Substituting (4.2) into (tO), we get 



^(z-2 + i)=e-2 + i. (4.3) 



This equation has two roots for ^: 



z2_2z + l 1 



6.2 = 1 + ± ^V(^-l)^[z2 + (4A2-2)z + l]. (4.4) 

Assume a boundary condition of the form 



fe=i 



10 



Substituting (4.2) into (4.5), we get 



1 1 " 

z ^ ^-^ 



To find Sfe, we have to find the Laurent expansion of the function on the left hand side. Let 



li{z) = V(z-l)2[z2 + (4A2-2)z + l]. 



(4.6) 



Observe that H{z) satisfies 



Hence 



H\z) = 2{z - \)\z^ + (3A2 - 2)^ + f - \^\lH{z). 



H'{z) ■ {{z ~ f )[z2 + (4A2 - 2)z + f]} = 2(2 - \){z'^ + (SA^ - 2)2 + f - \^\H{z). (4.7) 



Solving this equation by a Laurent series: tl(z) — fifnZ~"^, we obtain a recursion relation /z„ 

m>-2 

for m > 1, 

(771 + 2)^„ = [1 - 2X^ - m(4A2 - 3)]^„_i + [4 - GA^ - m(3 - 4A2)]//„_2 + (m - 3)^„„3, (4.8 



and 



/z_2 = 1, /i-i = 2A2 - 2, /io = 1 - 2A4 



(4.9) 



Then from (4.4) - (4.c), we have 



si = A^ 



s,^ = _l!±±^ for k>2. 



(4.10) 



(4.5) is nonlocal and has memory effects. In order to see how fast the memory decays, let us assume 



/ife ~ to" when to 2> 1, substituting into (4.5), and equating the coefficients of term order to", we get 



2 = (4A2 - 3)(1 + a) + (4 - GA^) + 2(3 - + a) + {2\^ - 2) - 3(1 + a). 



This gives a = — 2. The decay tendency of /ifc is shown in Figure 0. Here A = 0.01. 



If /C 7^ 0, we can proceed as before. But (4.4) changes to 
JC z^ ~2z+l 



1 



^.2 ±—^^[z^ + (/CA2 - 2)2 + 1][22 + (/CA2 + 4A2 - 2)z + 1]. (4.11) 
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Assuming a boundary condition of the form 



,n+l 



2ul - + A^K - (2 + /CX] + ^ 



(4.12) 



fe=0 



Substituting (4.2) into (4.12), we get 



1 1 

z-2 + --A2(--2-/C)==^ Skz~^. 
^ ^ k=i 

To find Sfe, we have to find the Laurent expansion of the function on the left hand side. Let g{lC, A) 
JCX^ + 2A2 - 2 and 



H{z) = v/[z2 + (/CA2 - 2)z + l][z2 + (/CA2 + 4A2 - 2)z + 1]. 



(4.13) 



Observe that H{z) satisfies 



H\z) = {2z^ + 3 giK., X)z^ + [g^{IC, A) + 2 - 4A'']z + .g(/C, X)}/H{z). 



Hence 



H'{z) • {z^ + 2 5(/C, A)z3 + [52(/^^ A) + 2 - AX'^]z^ + 2 g{K, \)z + 1} 
= i?(z)-{2z3 + 3g(/c,A)22 + [52(/C,A) + 2-4A-*]z + .g(/C,A)}. (4.14) 

Solving this equation by a Laurent series: H{z) — we obtain a recursion relation [irn 

m>-2 

for m > 2, 

(m + 2)^l,r^ - (2m + 1)[2 - \^{K: + 2)]^l„,^^ + (1 - m){2 - 4A'* + [2 - A2(/C + 2)f}ii„,^2 

+ (2 - /CA^ - 2A2)(2to - 5)Ai,„_3 + (4 - m)^l^-i, (4.15) 

and 

/x_2 = 1, M-i = ''CA^ + 2A2 - 2, Mo = 1 - aA'^, /ii = 2A'*(/CA2 + 2A2 - 2). (4.16) 



Then from ( |4ll| ) - ( jil^ ), we have 



So = —A /C, si = A , Sfe = — , for k>2 



(4.17) 
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In order to see how fast the memory decays, let us assume fj,k ^ m" when to ^ 1, substituting into 



(4.15), and equating the coefhcients of term order m" , we get 

= a{g{IC, A)(2 - a) - 2a[g^{IC, A) + 2 - AX^] - (6 + 9a)g(/C, A) - 8 - 8a}. 

This gives a = or a = -2/ |l + '^^^^+/^ A^}. 

These exact boundary conditions should be the same as the ones found numerically in Q]. It 



represents the exact Green's function for (3.1) which is nonlocal. However, this procedure appears to 
be impractical for realistic models, particularly when the atomistic region moves with time which is 
the case that interests us. But such calculations can at least give us guidelines on how to proceed to 
construct approximately reflectionless boundary conditions. 

4.2. Optimal Local Matching Conditions for ID Discrete Wave Equation 



A practical solution is to restrict (4J) to a finite number of terms and look for the coefficients 



{o-k,j} that minimize reflection. In order to do this, let us look for solutions of the type 



where lu is given by (|3.2D. R{^) is the reflection coefhcient at wavenumber ^. Inserting (4.18) into 



(4.1), we obtain 



The optimal coefficients {ofej } are obtained by 

min / WiC)\R{OfdC (4.20) 







subject to the constraint 

R{0) = 0,i?'(0) = 0, etc. (4.21) 

Here W{£,) is a weight function, which is chosen to be W{£,) = 1 in the examples below. 

Condition (4.21) guarantees that large scale information is transmitted accurately, whereas ( 4.2C| ) 



guarantees that the total amount of reflection is minimized. This procedure offers a lot of flexibility. 
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For example, instead of '^^ii minimize the total reflection over certain carefully selected 

interval. Another possibility is to choose the weight function to be the (empirically computed) energy 
spectrum. The coefficients {a^j} may then change in time to reflect the change of the nature of the 
small scales. In practice, we found it preferable to use * with some small 6, instead of 

Jo [^(OP'^'^j order to minimize the influence of ^ = tt for which we always have R{n) ~ 1. 



Let us look at a few examples. If in (4.1) we only keep the terms involving ai.o and ai^i, then 
imposing the condition i?(0) = gives 



iS = (1 - A^X"^ + Atul 



(4.22) 



If instead we keep terms involving ao,i, oi.o and oi^i, we can then impose both R{0) ~ and -R'(O) = 0. 
This gives us 

1 - 

(4.23) 



1 + At 



Conditions of the type (4.22) and (4.23) are intimately related to the absorbing boundary conditions 
proposed and analyzed in [^2| for the computation of waves. These conditions perform well for 
low wavenumbers but are less satisfactory at high wavenumbers. 

To improve the performance at high wavenumbers let us consider a case that include terms with 
k < 2, j < 3 and minimize ^ |i?(^)pc?^ (with 5 — 0.1257r) subject to the condition R{0) = 0, the 
optimal coefficients can be easily found numerically and are given by 



/ 



1.95264 



-7.4207 X 10" 



-0.95406 7.4904 x 10" 



-1.4903 X 10-2 ^ 



1.5621 X 10" 



(4.24) 



If instead we only include terms such that A: < 3, j < 2, then 

^ 2.9524 1.5150 x 10 



2 \ 



-2.9065 -3.0741 x 10" 



0.95406 1.5624 x 10-2 



(4.25) 



y V.HOWO i.0DZ4 X iU " J 

The resulting reflection coefficients R are displayed in Figure 
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4.3. Optimal Local Matching Conditions for Triangular Lattice 



The above procedure can be easily generalized. Let us take the triangular lattice as an example, 
and the boundary to be the x-axis. Given a boundary condition of the form: 



1<1 j 



(4.26) 



where Aj are some 2x2 matrices, and the summation is done over a pre-selected stencil, we can find 
the refiection matrix associated with this boundary condition. For that purpose, we look for solutions 
of the form 

"j" = E E Cpe'^^"^-"^-'^'^!]'!, (4.27) 
where a = 1,0 correspond to "incoming" and "outgoing" waves respectively (see Figure |^), /? = s,p 
correspond to "shear" and "pressure" waves respectively. Substituting into (^^) we obtain a relation 
between {C^,C^f and (Cf,C^)^, 



( r^O ^ 



Mo 



V J 



(cl\ 



= Mj 



(4.28) 



where 



Mo 

M/ 



-iuAt 



-iiuAt 



P?, C/p°]-5]^Aj [e^(«°"--J-'"^*)[/f, eK«°''--^-'"^*)[/0 

i<i j 



Zi^" Tj -lujAt) Jjl ^ ^iii'" Tj -luAt) Jjl 



In principle, we can solve the minimization problem 



min JwimMj^-MoiOfd^ 



(4.29) 
(4.30) 

(4.31) 



to find optimal {Aj'}, where the integration is over the Brillion zone. But in practice, we find it much 
more convenient to restrict the integration over a few selected low symmetry atomic planes. In the 
present context, it amounts to choosing special incidences where the phonons energy dominates. 
First, let us consider the case of normal incidence 9 = 90°. That means = 0. Then the matrix 



A in (3.4) becomes a diagonal matrix: 



n ^ 



2sin2^ 



15 



with two eigenvalues and eigenvectors: 



Ai = — sm — - — , Ui=l,0-', 
4 

108 . 2 Vi£.2a J. 

A2 = ^sin — - — , U2 = (0, 1)^ 



Then dispersion relations are 



WsAt = 2 arcsm( sm — — ), Us = (l,0) , 

ujp At — 2 a,rcsin( sm — - — ), Up = (0, 1) . 



(4.32) 
(4.33) 



If we take the absorb boundary condition as in ( 4.26 ), the matrices Mq and Mj are 



Mo 



Mr 



-ILuAt 



i<i j 



V 



(<i j 





/ 

(,i(&Vi-l'-^^t) 





(4.34) 



(4.35) 



where / is the 2x2 identity matrix. For consistency, we should require that the low wavenumber 



waves be transmitted accurately. Imposing (4.31), we get 



EEa], 



-lAt 







2a 



(4.36) 



(4.37) 



If we minimize (4.31) along normal incidence subject to the constraints (4.36) and (4.37), we obtain 



the desired matrices Aj'. For example, if we keep the terms with / = 0, 1 and j = (0,0), (—1,0), 



(—1, 1), the optimal coefficient matrices are 

/ 



^(0,0) 



ao - a" 

^(-1,0) - ^(-1,1) 



V 
/ 



0.947937634 -0.423061769£; - 09 
0.411523005^; - 09 0.911511476 

\ 



0.500000011 0.604865049£:- 08 
0.105260341E: - 07 0.499999996 
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■'^(-i.o) — -^(-1,1) 



V 



-0.473968784 -0.603638049i; - 08 
-0.102331855^; - 07 -0.455755718 
Next we consider the cases when both 6 = 60° and 6 = 120° are taken into account. For 6 = 60°, 
we have ^2 = V^^i, and 

9-4 sin^ ^ ^3(3 - 4 sin^ ^) 
^/3{3 - 4 sin^ ^) 15-12 sin^ ^ 
with two eigenvalues and eigenvectors: 

Ai = ^sin^ii^, Ui = (V3,-lf, 

A2 = f|sin2«i^(9-8sin2«i^), U2 = (l,\/3f. 
The dispersion relations are 



18 . 2^10 
A = sin — — 



2 arcsm sm — — 



ij„At = 2arcsin sin -^a/9 - 8sin^ 



Up = (1, Vsf. 



The consistency constraints are 
;<i j 



= Ai 



V 



V3 1 \ 
-1 V3 



EEa] 

1<1 j 



V3(C^" • rj /w - Z At ) Tj/LO-lAt 



for 61 = 60°, and 



I = EEa], 

1<1 j 



= Ai 



V 



EEa] 

1<1 j 



V3($^* • rj /w - / Ai) - {C^P Tj /w - « Ai) 
^ • rj/w - / Ai \/3(^^P • rj/w - I At) 



(4.38) 
(4.39) 



(4.40) 



(4.41) 



(4.42) 



(4.43) 



V3 -1 ^ 

1 Vs 

for 9 = 120°. For example, if we keep the terms for Z = 0, 1 and j = (0,0), (—1,0), (—1, 1), we have 
the optimal coefScient matrices 

/ 

0.929252841 -0.861918368E' - 09 

0.355336047^;- 09 0.908823412 



aO 

^(0,0) 



V 
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A(-1,0) 



^(-1,1) 



A(_i,o) 



V 
/ 

V 

f 

f 

V 



0.504255087 0.156041692£; - 01 
0.161793547£;- 01 0.499201553 



0.504255087 -0.156041692£; - 01 
0.161793547£; - 01 0.499201553 

\ 



-0.468881506 0.128308044£; - 01 
0.118075167^;- 01 -0.453613259 



-0.468881506 -0.128308044£; - 01 

-0.118075167£; - 01 -0.453613259 
If all three angles 9 = 60°, 90°, 120° are used with equal weight, then the optimal coefficient matrices 
are given by: 

/ 

0.963685659S+ 00 0.522045701£; - 05 
0.186532512S- 05 0.911580620£; + 00 



^(0,0) 



^(-1,0) 



ao - 

^(-1,1) - 



Ai 

^(-1,0) 



.A 

J 
\ 

J 



0.190155146£:+ 00 0.439544149£; - 02 
0.132862553^; - 01 0.497292487£; + 00 

0.190158427£:+00 -0.439289996£; - 02 
-0.132859770^;- 01 0.497292916^;+ 00 



-0.171943945£;+00 0.439598834£: - 02 
0.132858254^;- 01 -0.459575584£;+ 00 

-0.171944818£; + 00 -0.439293875£; - 02 
-0.132856732£; - 01 -0.453075576£; + 00 

5. ALGORITHMS AND IMPLEMENTATIONS 



The basic framework of our coupled continuum/atomistic method is that of an adaptive mesh 
refinement method p8[ | . The computational domain is covered by a grid that resolves the macroscopic 
features of problems, such as applied forces and boundary conditions. Regions near atomistic de- 
fects such as dislocations, interfaces, cracks, impurities, etc are detected using some error estimators. 
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Molecular dynamics are used in these regions to compute the location and momentum of each atom, 
together with the averaged quantities at the macroscopic grid points. Continuum equations are used 
elsewhere. At the interface between the two regions, matching conditions discussed in the last section 
are used. Specifically, we decompose the velocity and displacement fields into a large scale part and 
a small scale part. The large scale part is evolved using the values at the macroscopic grid points. In 
the atomistic regions, these are the averaged quantities. The small scale part is computed using the 
reflectionless boundary conditions discussed above. 

One important aspect of this method is the error estimators that are used to distinguish atomistic 
and continuum regions. The senstivity of the error estimators determines the balance between accuracy 
and efficiency. However, since there has already been a lot of work done on this specific problem 
| p9| , pO| pH , we will not pursue this question here further. We find it adequate in our work to use 
a refinement indicator (rather than an error estimator) which is given cither by an estimate of the 
stress, or a weighted average of the wavelet coefficients. 

Further details of our method are explained through a series of examples. 

5.1. Dislocation Dynamics in the Prenkel-Kontorova Model 

As the simplest model that encompasses most of the issues in a coupled atomistic/continuum 
simulation, we consider the Frenkel-Kontorova model 

Xj — Xj^i — 2xj + Xj-i — U'{xj) + f (5-1) 

where C/ is a periodic function with period 1, / is an external forcing. The continmmr limit of this 
equation is simply the Klein-Gordan equation 

utt = Uxx - Ku + f (5.2) 

where K — U"{0). We consider the case when there is a dislocation and study its dynamics under 
a constant applied forcing. We use U{x) = (x — [x])^ where [x] is the integer part of x. In this 



example we take (p3.1[) as our atomistic model, and (5.2) as our continuum model. For the coupled 



atomistic-continuum method, we use a standard second order finite difference method for (5^) in the 
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region away from the dislocation, and we use (5.1) in the region around the dislocation. However, we 



also place finite difference grid points in the atomistic region. At these points, the values are obtained 
through averaging the values from the atomistic model. At the interface between the atomistic and 
continuum regions, we decompose the displacement into a large scale and a small scale part. The large 
scale part is computed on the finite difference grid, using ( 4.25| ). The small scale part is computed 



using the reflectionless boundary conditions described earlier. The interfacial position between the 
MD and continuum regions is moved adaptively according to an analysis of the wavelet coefhcients 
or the local stress. The two strategies lead to similar results. Care has to be exercised in order to 
restrict the size of the atomistic region. For example, when wavelet coefficients are used in the criteria 
to move the atomistic region, we found it more efficient to use the intermediate levels of the wavelet 
coefhcients rather than the finest level. 

We first consider the case when a sharp transition is made between the atomistic and continuum 
regions with a f :f6 ratio for the size of the grids. Figure ^ is a comparison of the displacement and 
velocity fields computed using the full atomistic model and the coupled atomistic/continuum model, 
with / — 0.04. The atomistic region has 32 atoms. The full atomistic simulation has fO"^. Dislocation 
appears as a kink in the displacement field. Notice that at the atomistic/continuum interface, there is 
still substantial phonon energy which is then suppressed by the reflectionless boundary condition. No 
reflection of phonons back to the atomistic region is observed. In Figure ^, we compare the positions of 
the dislocation as a function of time, computed using the coupled method and the detailed molecular 
dynamics. Extremely good agreement is observed. 

We next consider a case with / — 0.02, which alone is too weak to move the dislocation, but to 
the left of the dislocation, we add a sinusoidal wave to the initial data. The dislocation moves as a 
consequence of the combined effect of the force and the interaction with the wave. Yet in this case 
the same atomistic/continuum method predicts an incorrect position for the dislocation, as shown in 
Figure |^. The discrepancy seems to grow slowly in time (see Figure p^ ). Improving the matching 
conditions does not seem to lead to significant improvement. 

The difference between this case and the case shown in Figure is that there is substantially more 

20 



energy at the intermediate scales. This is clearly shown in the energy spectrum that we computed 
for the two cases but it can also be seen in Figure ^ where an appreciable amount of small scale 
waves are present in front of the dislocation. Such intermediate scales are suppressed in a method 
that uses a sharp transition between the atomistic and continuum regions, unless we substantially 
increase the size of the atomistic region. We therefore consider the next alternative in which the 
atomistic/continuum transition is made gradually in a 1:2 or 1:4 ratio between neighboring grids. The 
right column in Figure || shows the results of such a method that uses a gradual 1:2 transition. We 
see that the correct dislocation position is now recovered. 

5.2. Friction between Flat and Rough Crystal Surfaces 

Our second example is the friction between crystal surfaces. To model this process atomistically, 
we use standard molecular dynamics with the Lennard- Jones potential |lj, |l^. First, we consider 
the case in which the two crystals are separated by a horizontal atomically flat interface. The atoms 
in the bottom crystal are assumed to be much heavier (by a factor of 10) than the atoms on top. To 
model the lack of chemical bonding between the atoms in the top and bottom crystals, the interaction 
forces are reduced by a factor of 5 between atoms in the top and bottom crystals. A constant shear 
stress is applied near the top surface. We use the periodic boundary condition in the x-direction. 

From a physical viewpoint, one interesting issue here is how dissipation takes place. Physically 
the kinetic energy of the small scales appears as phonons which then convert into heat and exit the 
system. A standard practice in modeling such a process is to add a friction term to the molecular 
dynamics in order to control the temperature of the system ||l4, 15 . In contrast, we ensure the proper 



dissipation of phonons to the environment by imposing the reflectionless boundary conditions for the 
phonons. The results presented below are computed using the last set of coefficient matrices presented 
at the end of Section 4. 

From Figure |ll] we see that we indeed obtain a linear relation between the mean displacement of 
the atoms in the top crystal as a function of time. The temperature of the system also saturates. 
Also plotted in Figure |l^ is the result of the mean displacement computed using the combined atom- 
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istic/continuum method. Here the continuum model is the Hnear elastic wave equation with Lame 
coefficients computed from the Lennard- Jones potential. The agreement between the full atomistic 
and the atomistic/continuum simulation is quite satisfactory. 

Next, we study the friction between two rough crystal surfaces. The setup is the same as before, 
except that the initial interface between the crystals takes the form y = f{x). The numerical results 
obtained are displayed in Figure In Figure we plot the positions of the atoms in the top and 
bottom crystals. We see that gaps are created in the case of rough interfaces. 

In Figure we compare the force- velocity relations for both flat and rough interfaces. Again the 
agreement between the coupled method and the full atomistic method is quite good. 

In the present problem, we used atomistic model in a narrow strip near the interface, and continuum 
model away from the interface. An interesting question is how wide the atomistic strip has to be. 
Clearly for the purpose of computational efficiency, we want the atomistic strip to be as narrow as 
possible. On the other hand, it has to be wide enough to provide an accurate description inside 
the boundary layer where important atomistic processes can be relaxed. There are two important 
atomistic processes in the present problem. The first is the vibration of the atoms around their local 
equilibrium positions. The second is the process of moving from one local equilibrium to the next, i.e. 
sliding by one atomic distance. Clearly the second process works on longer time scale. This process 
has to be resolved by the atomistic layer. In Figure ^5], we compare the atomic positions of a column 
of atoms which were initially vertical, i.e. they had the same a;-coordinates. From this picture one 
can also estimate the strain rate. We can clearly see that if the atomistic layer does not resolve the 
phonons generated by the second process, we get inaccurate results. 

5.3. Crack Propagation 



Our third example is the Slepyan model of fracture dynamics (2.11). In our coupled atom- 



istic/continuum method, we use full atomistic simulation (2.11) around the crack tip, and use ( ^.12| ) 



in the region far away from the crack tip. For the continuum equation, we use the displacement 
boundary condition u± — ±Un at the left boundary, and stress boundary condition ^ = at the 
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right boundary. Figure is a comparison of the fracture surface computed using the full atomistic 
model and the coupled atomistic/continuum method. 



Next we apply our method to the 2D Mode III fracture dynamics on a square lattice (2.14). Same 
boundary conditions as in the ID case are used for the continuum model. For the matching conditions 
between the atomistic and continuum regions, we used a stencil that consists of seven points: the values 
of the three nearest grid points next to the boundary at the current and previous time steps, plus 
the value at the boundary grid point at the previous time step. The optimization is carried out using 
angles = 45°, 90°, 135°. Figure ^ is a comparison of the fracture surface computed using the full 
atomistic model and the coupled atomistic/continuum method. Comparisons of the positions of the 
fracture tip as a function of time is given in Figure Esl The results are quite satisfactory. Finally in 



Figure 19, we display the shear waves generated as a result of the crack propagation, No reflection is 



seen. 



6. CONCLUSION 

In conclusion, we presented a new strategy for the matching conditions at the atomistic/continuum 
interface in multiscalc modeling of crystals. The main idea is to choose the boundary condition by 
minimizing the reflection of phonons along a few low symmetry atomic planes, subject to some accu- 
racy constraints at low wavenumbers. These conditions are adaptive if we choose the weight functions 



in (4.20) and (4.31) to reflect the evolving nature of the small scales. They minimize the reflection 
of phonons and at the same time ensure accurate passage of large scale information. The coupled 
atomistic/continuum method presented here is quite robust and works well at low temperature. At 
finite temperature and when nonlinearity is important at large scales, a new method has to be worked 
out. This work is in progress. 

We thank Tim Kaxiras for suggesting the problem of friction between rough interface. This work 
is supported in part by NSF through a PECASE award and by ONR grant N00014-01-1-0674. 
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List of Figures 



FIG. Triagular Lattice 

FIG. ^ 2D Slepyan model of fracture. The white dots indicate the equihbrium locations, the black 
dots indicate the displaced points once stress is applied. 

FIG. § Dispersion relation 

FIG. ^ Decay tendency of |^fc|. 



FIG. H Reflection coefficients for (4.24) and (4.25) 



FIG. Q The 'Incoming' and 'Outgoing' phonons near the boundary. 

FIG. Comparison of the displacement and velocity profiles computed using the full atomistic and 
the atomistic/continuum models, with / = 0.04. The top two graphs show the results in the whole 
computational domain. The bottom two graphs show the details near the dislocation. The solid line 
is the result of the atomistic/continuum method. The dash line is the result of the full atomistic 
method. 

FIG. § Comparison of the displacement and velocity profiles computed using the full atomistic and 
the atomistic/continuum models, with / = 0.02. The top two graphs show the results when the 
transition from the atomistic to continuum regions is sharp. The bottom two graphs show the results 
when the transition is gradual. Solid line is the result of the atomistic/continuum method. The dash 
line is the result of the full atomistic method. Only the region near the dislocation is shown. 

FIG. ^ Comparison of the positions of the dislocation as a function of time computed using the 
coupled method and the detailed molecular dynamics with / — 0.04. Dot line is the result of full MD 
simulation; solid line is the result with gradual transition between atomistic and continuum regions; 
dash line is the result with sharp transition. 
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FIG. 10 Comparison of the positions of the dislocation as a function of time computed using the 
coupled method and the detailed molecular dynamics with / = 0.02. Dot line is the result of full MD 
simulation; solid line is the result with gradual transition between atomistic and continuum regions; 
dash line is the result with sharp transition. 

FIG. Displacement and temperature as a function of time for the friction problem. 



FIG. 12 Displacement and temperature as a function of time for friction between rough surfaces. 



FIG. ^ The positions of the atoms near the interfaces. The white circles are light atoms, the black 
ones are heavy atoms. The top graph is the initial state, the bottom graph is the late state at t = 1000. 



FIG. 14 Comparison of the force-velocity relations for both flat and rough interfaces. The top two 
lines are the results for flat case, the bottom two lines are the results for rough case. The solid lines 
are the results for coupled atomistic/continuum method, the dash lines are the results for full MD 
simulation. 



FIG. 15 Comparison of the atomic positions of a column of atoms which had the same a;-coordinates 
initially. Solid line is the result of full MD, the line with o is the result of coupled method with 96 
layers in the atomistic region, the line with + is the result of coupled method with 16 layers in the 
atomistic region. In the coupled method, the ratio of atomistic and continuum grids is 1:8. 

FIG. ^ One-dimensional fracture problem. The left graph shows the fracture surface at time t — 
and the right one shows the fracture surface at time t = 600 with b — 0.01, N — Un = 4. The dash 
line is the result of the full molecular dynamics simulation. The solid line is the result of the coupled 
atomistic/continuum method. The ratio of atomistic and continuum grids is 1:8. 

FIG. ^ Two-dimensional fracture problem. The left graph shows the fracture surface at time t = 
and the right one shows the fracture surface at time t = 200 with b = 0.01, = 512, and Un = VN- 
There are 800 atoms in each row. The dash line is the result of the full molecular dynamics simulation. 
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The solid line is the result of the coupled atomistic/continuum method. In coupled method, we divide 
the whole domain into three parts. The middle part including the crack surface with 800 x 64 atoms 
is the MD region. The top and bottom parts are continuum regions. The ratio of atomistic and 
continuum grids in each dimension is 1:8. 



FIG. 18 Comparisons of the positions of the crack-tip as a function of time. The dash line is the 



result of the full MD simulation, the solid line is the result of the coupled method. 



FIG. 19 The shear waves. We divide the whole domain with 2048 x 2048 atoms into three parts. 
The middle part including the crack surface with 2048 x 128 atoms is the MD region. The top and 
bottom parts are continuum regions with 128 x 62 finite difference grids in each region. 
FIG. Efl The enlarged picture of the MD region near the crack-tip. 
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